#include <feel/options.hpp>
#include <feel/feelcore/application.hpp>

#include <feel/feelalg/backend.hpp>

#include <feel/feeldiscr/functionspace.hpp>

//#include <feel/feelpoly/im.hpp>

#include <feel/feelfilters/gmsh.hpp>

#include <feel/feelfilters/exporter.hpp>
//#include <feel/feelpoly/lagrange.hpp>
//#include <feel/feelpoly/crouzeixraviart.hpp>

//#include <feel/feelmesh/elements.hpp>

#include <feel/feelvf/vf.hpp>

#include <feel/feelfilters/geotool.hpp>
/**
 * This routine returns the list of options using the
 * boost::program_options library. The data returned is typically used
 * as an argument of a Feel::Application subclass.
 *
 * \return the list of options
 */

namespace Feel
{
namespace Cavity
{

inline
Feel::po::options_description
makeOptions()
{
    Feel::po::options_description stokesoptions("Stokes options");
    stokesoptions.add_options()
        ("penal", Feel::po::value<double>()->default_value( 0.5 ), "penalisation parameter")
        ("f", Feel::po::value<double>()->default_value( 0 ), "forcing term")
        ("mu", Feel::po::value<double>()->default_value( 1.0 ), "reaction coefficient component")
        ("hsize", Feel::po::value<double>()->default_value( 0.1 ), "first h value to start convergence")
        ("bctype", Feel::po::value<int>()->default_value( 0 ), "0 = strong Dirichlet, 1 = weak Dirichlet")
        ("bccoeff", Feel::po::value<double>()->default_value( 100.0 ), "coeff for weak Dirichlet conditions")
        ("export-matlab", "export matrix and vectors in matlab" )
        ;
    return stokesoptions.add( Feel::feel_options() ) ;
}


/**
 * This routine defines some information about the application like
 * authors, version, or name of the application. The data returned is
 * typically used as an argument of a Feel::Application subclass.
 *
 * \return some data about the application.
 */
inline
Feel::AboutData
makeAbout()
{
    Feel::AboutData about( "stokes_cavity" ,
                           "stokes_cavity" ,
                           "0.1",
                           "Stokes equation on simplices or simplex products",
                           Feel::AboutData::License_GPL,
                           "Copyright (c) 2009 Universite de Grenoble 1 (Joseph Fourier)");

    about.addAuthor("Christophe Prud'homme", "developer", "christophe.prudhomme@ujf-grenoble.fr", "");
   return about;

}


/**
 * \class Stokes class
 * \brief solves the stokes equations
 *
 */
template<int Dim,
         typename BasisU,
         typename BasisP,
         template<uint16_type,uint16_type,uint16_type> class Entity>
class Stokes
    :
        public Application
{
    typedef Application super;
public:


    typedef double value_type;

    typedef Backend<value_type> backend_type;
    typedef boost::shared_ptr<backend_type> backend_ptrtype;

    /*matrix*/
    typedef typename backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
    typedef typename backend_type::vector_ptrtype vector_ptrtype;

    /*mesh*/
    typedef Entity<Dim,1,Dim> convex_type;
    typedef Mesh<convex_type> mesh_type;
    typedef boost::shared_ptr<mesh_type> mesh_ptrtype;

    /*basis*/
    typedef BasisU basis_u_type;
    typedef BasisP basis_p_type;
    typedef Lagrange<0, Scalar> basis_l_type;
    typedef bases<basis_u_type,basis_p_type, basis_l_type> basis_type;

    /*space*/
    typedef FunctionSpace<mesh_type, basis_type> space_type;
    typedef boost::shared_ptr<space_type> space_ptrtype;

    BOOST_MPL_ASSERT( ( boost::is_same<typename space_type::bases_list, basis_type> ) );
    BOOST_MPL_ASSERT( ( boost::is_same<typename mpl::at<typename space_type::bases_list,mpl::int_<0> >::type, basis_u_type> ) );
    BOOST_MPL_ASSERT( ( boost::is_same<typename mpl::at<typename space_type::bases_list,mpl::int_<1> >::type, basis_p_type> ) );
    BOOST_MPL_ASSERT( ( boost::is_same<typename mpl::at<typename space_type::bases_list,mpl::int_<2> >::type, basis_l_type> ) );

    /* functions */
    typedef typename space_type::element_type element_type;
    typedef typename element_type::template sub_element<0>::type element_0_type;
    typedef typename element_type::template sub_element<1>::type element_1_type;
    typedef typename element_type::template sub_element<2>::type element_2_type;

    /* export */
    typedef Exporter<mesh_type> export_type;
    typedef boost::shared_ptr<export_type> export_ptrtype;

    Stokes( int argc, char** argv, AboutData const& ad, po::options_description const& od )
        :
        super( argc, argv, ad, od ),
        M_backend( backend_type::build( this->vm() ) ),
        meshSize( this->vm()["hsize"].template as<double>() ),
        exporter( Exporter<mesh_type>::New( this->vm(), this->about().appName() ) )
    {
        mu = this->vm()["mu"].template as<value_type>();
        penalbc = this->vm()["bccoeff"].template as<value_type>();
    }


    /**
     * create the mesh using mesh size \c meshSize
     */
    mesh_ptrtype createMesh( double meshSize );


    /**
     * run the convergence test
     */
    void run();

private:

    /**
     * solve system
     */
    void solve( sparse_matrix_ptrtype const& D, element_type& u, vector_ptrtype const& F, bool is_sym );

    /**
     * export results to ensight format (enabled by  --export cmd line options)
     */
    void exportResults( element_type& u );

private:

    backend_ptrtype M_backend;
    double meshSize;

    double mu;
    double penalbc;

    boost::shared_ptr<export_type> exporter;
}; // Stokes

template<int Dim, typename BasisU, typename BasisP, template<uint16_type,uint16_type,uint16_type> class Entity>
typename Stokes<Dim,BasisU,BasisP,Entity>::mesh_ptrtype
Stokes<Dim,BasisU,BasisP,Entity>::createMesh( double meshSize )
{
    mesh_ptrtype mesh( new mesh_type );

    GeoTool::Node x1(0,0);
    GeoTool::Node x2(1,1);
    GeoTool::Rectangle R( meshSize,"LEFT",x1,x2);
    R.setMarker(_type="line",_name="Lift",_marker3=true);
    R.setMarker(_type="line",_name="Paroi",_marker1=true,_marker2=true,_marker4=true);
    R.setMarker(_type="surface",_name="OmegaFluide",_marker1=true);

    GeoTool::Node x_centre(0.5,0.5);
    GeoTool::Node x_bord(0.5,0.7);
    GeoTool::Circle C( meshSize/6.,"disque",x_centre,x_bord);
    C.setMarker(_type="line",_name="Paroi",_markerAll=true);
    C.setMarker(_type="surface",_name="OmegaDisque",_marker1=true);

#if 0
    //GmshTensorizedDomain<convex_type::nDim,convex_type::nOrder,convex_type::nRealDim,Entity> td;
    //td.setCharacteristicLength( meshSize );
    //std::string fname = td.generate( convex_type::name().c_str() );
    Gmsh gmsh;
    std::string fname = gmsh.generate( "domain", (R-C).geoStr(),false,false,false  );

    ImporterGmsh<mesh_type> import( fname );
    mesh->accept( import );
    return mesh;
#endif
    return (R-C).createMesh<mesh_type>("domain");
} // Stokes::createMesh


template<int Dim, typename BasisU, typename BasisP, template<uint16_type,uint16_type,uint16_type> class Entity>
void
Stokes<Dim, BasisU, BasisP, Entity>::run()
{
    if ( this->vm().count( "help" ) )
        {
            std::cout << this->optionsDescription() << "\n";
            return;
        }

    using namespace Feel::vf;

    this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                            % this->about().appName()
                            % convex_type::name()
                            % BasisU::nOrder
                            % this->vm()["hsize"].template as<double>()
                            );

    mesh_ptrtype mesh = createMesh( meshSize );

    space_ptrtype Xh = space_type::New( mesh );

    element_type U( Xh, "u" );
    element_type V( Xh, "v" );
    element_0_type u = U.template element<0>();
    element_0_type v = V.template element<0>();
    element_1_type p = U.template element<1>();
    element_1_type q = V.template element<1>();
    element_2_type lambda = U.template element<2>();
    element_2_type nu = V.template element<2>();

    vector_ptrtype F( M_backend->newVector( Xh ) );

    AUTO( deft, gradt(u) );
    AUTO( def, grad(v) );

    // total stress tensor (trial)
    AUTO( SigmaNt, (-idt(p)*N()+mu*deft*N()) );

    // total stress tensor (test)
    AUTO( SigmaN, (-id(p)*N()+mu*def*N()) );

    // u exact solution
    AUTO( u_exact, vec(cos(Px())*cos(Py()),sin(Px())*sin(Py()) ) );

    // this is the exact solution which has zero mean : the mean of
    // cos(x)*sin(y) is sin(1)*(1-cos(1))) on [0,1]^2
    AUTO( p_exact, cos(Px())*sin(Py())-(sin(1.0)*(1.-cos(1.0)) ) );

    // f is such that f = \Delta u_exact + \nabla p_exact
    //AUTO( f, vec(2*cos(Px())*cos(Py())-sin(Px())*sin(Py()),
    //             2*sin(Px())*sin(Py())+cos(Px())*cos(Py()) ) );

    auto f=vec(cst(0.),cst(0.));

    auto vit=vec( Px()*-(1-Px()) , cst(0.));

    // right hand side
    form1( Xh, F, _init=true )  =
        integrate( elements(mesh), trans(f)*id(v) )+
        integrate( markedfaces(mesh,mesh->markerName("Lift")), trans(vit)*(-SigmaN+penalbc*id(v)/hFace() ) );

    Log() << "[stokes] vector local assembly done\n";

    sparse_matrix_ptrtype D( M_backend->newMatrix( Xh, Xh ) );

    form2( Xh, Xh, D, _init=true )=integrate( elements(mesh), mu*trace(deft*trans(def)) );
    form2( Xh, Xh, D )+=integrate( elements(mesh), - div(v)*idt(p) + divt(u)*id(q) );
    form2( Xh, Xh, D )+=integrate( elements(mesh), id(q)*idt(lambda) + idt(p)*id(nu) );
    form2( Xh, Xh, D )+=integrate( boundaryfaces(mesh),
                                   -trans(SigmaNt)*id(v)
                                   -trans(SigmaN)*idt(u)
                                   +penalbc*trans(idt(u))*id(v)/hFace() );

    D->close();
    F->close();

    this->solve( D, U, F, false );

    //v = vf::project( Xh->template functionSpace<0>(), elements( Xh->mesh() ), u_exact );
    //q = vf::project( Xh->template functionSpace<1>(), elements( Xh->mesh() ), p_exact );

    this->exportResults( U);

} // Stokes::run

template<int Dim, typename BasisU, typename BasisP, template<uint16_type,uint16_type,uint16_type> class Entity>
void
Stokes<Dim, BasisU, BasisP, Entity>::solve( sparse_matrix_ptrtype const& D,
                                   element_type& u,
                                   vector_ptrtype const& F,
                                   bool is_sym )
{
    vector_ptrtype U( M_backend->newVector( u.functionSpace() ) );
    M_backend->solve( D, D, U, F, false );
    u = *U;
} // Stokes::solve

template<int Dim, typename BasisU, typename BasisP, template<uint16_type,uint16_type,uint16_type> class Entity>
void
Stokes<Dim, BasisU, BasisP, Entity>::exportResults( element_type& U )
{

    export_ptrtype __exporter( export_type::New( "ensight",//"gmsh",
                                                 this->about().appName() ) );
                               //+ "_"
                               //                + mesh_type::shape_type::name() ) );

    // if ( exporter->doExport() )
    // {
    //__exporter->step( 0 )->setMesh( U.template element<0>().functionSpace()->mesh() );
        __exporter->step( 0 )->setMesh( U.functionSpace()->mesh() );
        __exporter->step( 0 )->add( "u", U.template element<0>() );

     __exporter->step( 0 )->add( "p", U.template element<1>() );
        //        exporter->step( 0 )->add( "u_exact", V.template element<0>() );
        //exporter->step( 0 )->add( "p_exact", V.template element<1>() );
    __exporter->save();
    //}
} // Stokes::export
} // Feel
} // Cavity
int
main( int argc, char** argv )
{

    using namespace Feel;
    /* assertions handling */
    Feel::Assert::setLog( "stokes.assert");

    const int nDim = 2;
    typedef Cavity::Stokes<nDim, Lagrange<2, Vectorial>,Lagrange<1, Scalar>, Simplex> stokes_type;
    //typedef Feel::Stokes<nDim, CrouzeixRaviart<1, Vectorial>,Lagrange<0, Scalar>, Simplex> stokes_type;


    /* define and run application */
    stokes_type stokes( argc, argv, Cavity::makeAbout(), Cavity::makeOptions() );
    stokes.run();
}





